

forv yearval=2001/2010{
*****************
* CENSORED MLE *
*****************
forv eta=1/6{
* manipulation window 
scalar lftcff1=0.15-0.005*(ceil(`eta'/3)+1)   
scalar rtcff1=0.15+0.005*(ceil(`eta'/3)+1)*(1+`eta'-floor((`eta'-1)/3)*3) 
use "./data/pm10data.dta", clear
qui keep if id==idval // city by city
qui keep if year==`yearval' // year by year
scalar cff1=0.15
sum pmconc if pmconc<=float(cff1)
gen Fxnp=r(N)/_N  
gen cens=0
replace cens=1 if pmconc>float(lftcff1) & pmconc <=float(cff1)  
gb2lfit_mod pmconc, cdf(pmparacdf1)  censvar(cens)
local ahat=e(ba)
local bhat=e(bb)
local phat=e(bp)
local qhat=e(bq)
local atx=cff1
gen Fxp=ibeta(`phat',`qhat', (`atx'/`bhat')^`ahat'/(1+(`atx'/`bhat')^`ahat'))
qui gen propFx=(Fxnp-Fxp)/Fxnp*100

keep city year id propFx
gen rtcff1=rtcff1
gen lftcff1=lftcff1
keep in 1
local cityname=city[1] 

if `yearval'==2001 & `eta'==1{
save "./data/`cityname'_windowsensitivity.dta", replace
}
else {
append using "./data/`cityname'_windowsensitivity.dta"
save "./data/`cityname'_windowsensitivity.dta", replace
}
}

}


twoway (connected  propFx year if lftcff1==float(0.135) & rtcff1==float(0.18))  ///
(line  propFx year if lftcff1==float(0.14) & rtcff1==float(0.17), lcolor(gs1) lpattern(dot)) ///
(line propFx year if lftcff1==float(0.14) & rtcff1==float(0.18), lcolor(gs1) lpattern(dot)) ///
(line propFx year if lftcff1==float(0.14) & rtcff1==float(0.19), lcolor(gs1) lpattern(dot)) ///
(line propFx year if lftcff1==float(0.135) & rtcff1==float(0.195), lcolor(gs1) lpattern(dot)) ///
 (line propFx year if lftcff1==float(0.135) & rtcff1==float(0.21), lcolor(gs1) lpattern(dot)), ///
ytitle("") legend(order(1 "Benchmark" 2 "Robustness Checks")) ///
xlabel(2002(2)2010) ylabel(`foo') graphregion(color(white)) ///
name(robust, replace) nodraw title("Robustness Check")




